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Abstract. We develop an analytic and environment-dependent interatomic potential for the 
overlap repulsion in solid argon, based on an approximate treatment of the non-orthogonal 
Tight-Binding theory for the closed-shell systems. The present model can well reproduce 
the observed elastic properties of solid argon including Cauchy violation at high pressures, 
yet very simple. A useful and novel analysis is given to show how the elastic properties are 
related to the environment-dependence incorporated into a generic pairwise potential. The 
present study has a close link to the broad field of computational materials science, in which 
the inclusion of environment dependence in short-ranged repulsive part of a potential model is 
sometimes crucial in predicting the elastic properties correctly. 

PACS numbers: 62.50.+p, 62.20.Dc 

1. Introduction 

Recent progress of Brillouin spectroscopy at very high pressuresflU 13 has revealed that 
interatomic forces in fee solid argon must be far beyond any kinds of two-body, central 
force model. Shimizu et al. flU precisely measured a large violation of the Cauchy relation 
up to 70 GPa and stressed the important role of many-body forces, in order to construct 
good potentials for high-density noble gases, which should be crucial in understanding their 
behavior in planetary bodies by means of molecular simulations. 

The Cauchy relation[[3|| for the elastic constants of cubic crystals at a hydrostatic pressure 
P is given by C\ 2 — C44 — 2P = 0, which must be satisfied in centrosymmetric cubic 
crystals, including the fee solid argon, if the total energy is given by sum of purely pairwise 
terms. The deviation from it is therefore a measure of the many-atom nature of interatomic 
interactions. The Aziz-Slaman model for high pressure argon[4], which might be one of 
the most sophisticated yet simple ones thus far proposed, fails to reproduce any violation of 
Cauchy relation, simply because the model is pairwise. 

On the other hand, the ab initio Density Functional Theory (DFT) approach 
using the pseudopotential planewave method|[5]|, and that with projector- augmented wave 
implementation for core electrons [6], and linear muffin-tin orbital (LMTO) methodEUHl have 
successfully reproduced the observed elastic constants, as well as the density of the solid argon 
over the measured range of the pressure. From these theoretical results, it should be a natural 
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Figure 1. Total Cauchy violation S, decomposed into contributions from the kinetic (thin 
solid), electrostatic (dot-dashed), and exchange-correlation (dashed) energies. 



and logical consequence that one would expect to have an simple and reasonably accurate 
model for the many-atom interaction in condensed argon, by coarse-graining the ab initio 
electronic models to a rather empirical atomistic model. 

The importance of many-body forces in solid argon at high pressures has been examined 
by several authors along the idea of many-atom expansion[[9l[T0l[TT||, m which one assumes 
that the total energy is 'additively' decomposable into iV-atom (N = 2, 3, 4 • ■ ■) terms plus 
the additional energy of zero-point vibrations, and that the expansion is well convergent. The 
three-atom contribution from the exchange energy was emphasized[10] because it stabilise 
argon in fee structure, rather than the hep, which is predicted by all the available pairwise 
models without the zero-point energy [fTTT] . However, it is pointed out that the convergence 
of this type of expansion becomes worse in a situation in which many-atom effect is more 
import ant |fl"Tfl. 

Figured] shows the Cauchy violation defined by 

5 = C 12 - C 44 - 2P, (1) 

and its breakdown into the contributions from kinetic, electrostatic and exchange-correlation 
energies predictedflU by using all-electron calculation within DFT lfT2l . Each curve is plotted 
as a function of the total pressure. Clearly, the central role for the observed negative 5 is played 
by kinetic energy, which remains after the large cancellation by the opposite contributions 
from electrostatic and exchange-correlation. It should be noted that non-zero contribution 
from the electrostatic energy immediately excludes a primitive picture of overlapping frozen 
atomic charge-density. Therefore, it is implied that deformation of density (wavefunctions) 
should be relevant. 

The purpose of this paper is to develop an analytic interatomic potential for solid argon, 
that is based on the quantum mechanics of electrons, and that can well reproduce the observed 
elastic properties including the Cauchy violation at high pressures, yet that is made as 
simple as possible. The problem we are to treat now has a close link to the broad field of 
computational materials science, since the inclusion of environment dependence in short- 
ranged repulsion is sometimes crucial |[i~3l [i~4l [i~5l [i~6ll in obtaining reliable and transferable 
models for simulations in empirical and semi-empirical approaches. 



Environment-dependent potential for solid Ar 



3 



In section 2, a Tight-Binding description for overlap repulsion in closed-shell atoms, 
which depends on atomic environment, will be presented as the theoretical base that 
underpins more empirical and analytic model. General properties of repulsive potential with 
environment-dependent parameters are analysed, and a simple functional form for the overlap 
repulsion is designed for argon and proposed in section 3. The result of the fitted analytic 
potential is presented in section 4, followed by a concluding section 5. 



2. Environment-dependent overlap repulsion: Tight-Binding description 

A system of closed-shell atoms may be well described within the non-orthogonal Tight- 
Binding Bond Model (TBBM) ifTTl [T8l . in which the total binding energy is given by 

Eb = E cov + E ren + E mp + -Evdw- (2) 
The first term is the covalent energy 

E cov = 2Tr[HS- 1 ] - 2Tr[H] = -^[HOS -1 ], (3) 
where H and S = 1 + are the Hamiltonian and overlap matrices defined by 



(H) Wv = / C(r)ffiMr)d 8 r W 



and 



(S)i M v = ^m<W + (°)ifi,ju = J ^>( r )^v(r)d 3 r, (5) 

in a basis set of atomic orbitals Tp ifJr , where p runs over orbitals on site i. The spin degeneracy 
enters as the prefactor 2 before the usual symbol Tr for the trace. The Hamiltonian operator H 
refers to the density of superposition of frozen atomic densities [fT8l IT7T1 . Note that the inverse 
overlap matrix, S -1 , in Eq.© is equivalent to the density matrix in this case of fully occupied 
system. The second term in Eq. ©, E ven , is defined by 

£ ren = 2Tr[H]-2Tr[H fa ], (6) 

which accounts for the on-site energy shift due to the contraction or localisation of free atomic 
orbitals on going into a condensed environment, and H fa is a diagonal matrix of the free 
atomic energy levels. E Tep represents the contribution from the change in the electrostatic 
and exchange-correlation energies associated with frozen atomic charge densities as they are 
brought together from the free space. Thus, this term is environmentally independent by 
construction, and usually approximated as a sum of repulsive pair-wise potentials between 
atoms. E vd y?, added supplementary to TBBM, denotes the van der Waals potential, which 
may be approximated using the empirical pair-wise inverse power function of interatomic 
separation, that is — c 6 i?~ 6 . 

In order to illustrate that Eq.© essentially represents the overlap repulsion, we now 
consider only the outermost p-shell states as the basis, and then the bond integral[19J part, B, 
may be separated from H as 

H = e p S + B, (7) 
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where e p is the common diagonal element of H. Further simplification can be made exploiting 
the Wolf sberg-Helmholtz |fT9ll (or extended Hiickel) approximation to write B = —bO using a 
constant b > 0. We easily find 

£ cov = 2Wr[0 2 (l + 0)- 1 ], (8) 

and see immediately that E cov = if the basis is orthogonal (i.e. = zero), and it turns into 
repulsive when the matrix is switched on. Note that Eq.® has no explicit dependence on 
e p . The name of this term, therefore, is only nominal for the noble gases as it gives overlap 
repulsion rather than the covalent bonding. The lowest order term in Eq. ([8]) is given by 
26Tr [O 2 ], which may be broken down into contributions in a purely pair- wise form 

= Ab{\O ppa (R tJ )\ 2 + 2|O pp7r (%)| 2 }, (9) 

where O ppo -(Rij) and O pp7r (Rij) are a- and 7r-overlap integrals along interatomic separation 
Hij. The higher order correction terms, which arise from multiplication of (l+O) -1 in Eq. ([8]), 
or alternatively, multiplication of (1 + O) -1 / 2 from both sides of O 2 in a symmetric Lowdin 
form, would introduce many-atom effects as derived by Nguyen-Manh et al. for environment- 
dependent bond integrals [20]. For simplicity of our model, we may neglect these higher order 
corrections to write 

E cov Tr[-2HO] = \ ]T $y (i^) = E ovh (10) 

Assuming the Slater-type atomic p-orbitals with exponential tail of exp(— k^) for atom i, 
the bond and overlap integrals decay like exp[—fa + Kj)Rij] and the overlap potential Eq. © 
takes form of 

^ij(Rij) = (polynomial of R^) x exp[— fa + Kj)R\j\ (11) 

An important environment effect can naturally be taken into account if we think of fa} as a 
set of variational parameters. It is a straightforward exercise to show for a hydrogen-like atom 
with an effective atomic number Z* that (in atomic Rydberg units) 

J W(r)(-V 2 - 2Z7r)W(r)d 3 r = fa - Z* /2) 2 - (Z* /2) 2 (12) 

with —(Z*/2) 2 being the lowest p-energy level, we see that the parabolic penalty for an 
augmented k arises from increase in kinetic energy due to localisation. This parabolic 
behaviour occurs as a result of the change in the effective radius of the atomic wavefunction, 
regardless of particular form of the atomic pseudopotentials. Therefore, we may assume for 
the diagonal matrix elements of H that 

(H)i Mi i M = fa — k m ) 2 + e i0 = e ip fa), (13) 

where n i0 and e i0 are constants, and thus, E ien has the role of penalty for localisation of 
atomic orbitals through Eq. (fT3l . while the localisation will reduce the magnitude of overlap 
repulsion. Therefore the optimum values of will be determined by minimizing E ov i+E Ten 
with respect to each Ki. In the case of p- shells, given a simplified form for overlap potential 
$ij = 6gexp[— fa + Kj)Rij) with a constant q, this minimisation leads to 

AKi = Ki- K i0 = 1 R ij ex P[-( K i + 04) 

i(^) 
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The set of solutions is indeed environment-dependent and we see that the constant K i0 
is the solution for the limiting case of infinitely separated atoms. Since the penalty due to 
increase in the kinetic energy is steep, AKi/n i0 would be small enough to replace /% in the 
exponential in Eq.(fT4l) with K i0 . This solves the equation to give explicit A«j that is exact to 
first order. The energy increase AE rcn due to this minimisation is given exactly by the sum of 
3(AKj) 2 , which partially sets off and just halves the lowest order decrease in overlap energy 
AE ov \. The resultant lowering, ^AE ovh could be obtained by employing |A/«j instead of A/% 
in the overlap energy. This treatment eliminates AE rcn and simplifies the functional form of 
the potential, which is to be proposed in the next section. 

The environmental effect that we are looking at is a tendency that more contracted atomic 
orbitals are preferred in a denser environment. The physics behind it has been beautifully 
justified in the pioneering work by Skinner and Pettifor[|2T|. who have implemented the 
chemical pseudopotential theory using the orbital exponent as a variational parameter within 
the Harris-Foulkes scheme Il22ll23l . and found that the orbital exponents (k's in our notation) 
of hydrogen atoms in molecule, simple cubic and fee lattices are strongly environment- 
dependent as stated above. 



3. Analytic model for solid argon 

Let us first analyse some general properties of a repulsive potential that is a function of 
environment-dependent parameters as well as the distance. 

Provided that the functional form of repulsive interatomic potential is given by 

Qi^RijlXi + Xj) (15) 

with the environment-dependent parameters (A's) that are written as a sum of pairwise 
functions, namely 

Xi = E Pi R ik)- (16) 

fc(#) 

Cauchy violation can be calculated analytically for cubic lattice. The result is writtenjHl 



2 

9ft 



(17) 



a = E Rkp'(Rk), (18) 

k&O) 



with 



where Q is the atomic volume and Rj = i? j — y + y] + z] is the distance to atom j from 
the central atom i = at the origin. The prime on p denotes the distance derivative. Note 
that all lattice sites are equivalent under a homogeneous strain. a represents the strength 
of environmental effect on atom 0. Since p(R) at this point is completely arbitrary, we may 
assume that it is a positive and monotonically decreasing function of distance in the range of 
interest; hence a < 0. It may also be a physically reasonable assumption that the repulsive 
potential $oj(> 0) is also a monotonically decreasing function of distance in the range of 
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interest. We see from Eq.(fT71) for the case of very weak a that negative Cauchy violation 
occurs if d 2 $ j/dRjd\ > 0. This condition is likely to be fulfilled, since an environmental 
effect tends to weaken the overlap repulsion to give d&oj/dXo < 0, as we have seen in the 
previous section. The expression for the pressure from $ is given by 



p m 



(19) 



The first term in the square bracket represents the pairwise component. We see environmental 
effect, the second term, causes reduction in pressure as expected. 

A simple functional form of overlap repulsion that takes into account the physics of the 
environment effect as we have discussed is now proposed, that is, 

fcyCRy; Xi + Xj) = exp(-Ai) exp(-A i )V R ( J R ij ), (20) 

where Vr is environmentally independent pairwise function. The environmental effect on site 
i is entering in a very simple separable form by a factor exp(— A»), which corresponds to the 
contraction factor exp(— AniRij). However, the direct dependence on the particular length 
Rij has been dropped for simplicity. This manner of parameterisation for the environmental 
effect can also be seen in the 'breathing-shell model' (See Ref. E4l |25l and references 
therein) and 'compressible ion model' ll26ll for oxides, such as MgO. The both models are 
provided with the parameters that correspond to the variation in effective size of ionic cores, 
and reduction in core size causes exponential reduction as Aj in the present model does. It 
should be noted that a kind of penalty function has been eliminated from the present model 
as it was justified in the previous section. The contraction factor of type exp(— A«ji2y) was 
modelled by Nguyen-Manh et al. in the form of screened Yukawa-type potential [fT6l[T5l and 
it was used to explain Cauchy pressures in transition metals intermetallics within a Tight- 
Binding and Harris-Foulkes approaches. 

It follows from the functional form proposed above that full expressions for the pressure 
P, Cauchy violation 5, adiabatic bulk modulus B, and the cubic elastic constants Cn, Ci 2 , C 44 
are given by 

P =i(-u + 2ua ), (21) 
5 =^(-a v + ual), (22) 
B =lp + lK + 5, (23) 
C u = - P + P s + K s + 6, (24) 
C12 = |(3P + K - P s - K s ) + 5, (25) 
C 44 = \{-P + K - P s - K s ) (26) 



with 



P s = \(-v s + 2ua s Q ), (27) 



K = \{w - 2u/3 ), K s = \{w s - 2u(3 s ), (28) 



3l v -r^«<"07> 
\{w - 2up ), h - 3V 

where u equals the energy density and v, w also are similar quantities determined by first- and 
second-order derivatives of the potential, and u s ,v s , w s are weighted sums. These are defined 
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Figure 2. The Cauchy violation, 6, and pressure, P, as a function of strength of environment- 
dependence, cto, for the repulsive potential. 



by 

u = 4t E ^ « S = ^ E $ oi*i, (29) 
«=iE ^o;, " 8 = ^E ^0^, (30) 

™ = ^£*X-> - S = ^E^ 2 ^ (3D 

j(#o) j(^o) 

with Sj = (xj + yj + zj)/Rlj. Together with a , three of other quantities representing the 
strength of environment-dependence are defined: 

a s = J2 Rkp'(Rk)s k , (32) 

fc(^0) 

A) = E Rlp"(Rk), Po = E Rlp"(Rk)sk- (33) 

fc(/0) fc(^0) 

Equations ((211) and (|22l) . seen as a linear and quadratic functions of a respectively, 
explain how the negative Cauchy violation occurs when an environment-dependence is 
introduced, as presented instructively in Fig.[2l It is to be noted that ceo = v/ (2w) is unphysical 
point where our 'repulsive' potential is found no longer repulsive, giving P = 0. Therefore, 
the magnitude of dimensionless parameter a should usually be much smaller than \v/(2u)\. 
An instructive example may be the case of inverse-power-law potential, $ oc R~ n , in which 
the critical value can be easily found to be v/(2u) = —n/2. These analysis should be useful 
in understanding how the environmental dependence in repulsive worked for the problem of 
small or negative Cauchy pressure (Ci 2 — C 44 for cubic systems, C ri — C 44 and Ci 2 — C 66 for 
tetragonal or hexagonal systems) in transition metals and intermetallic compounds lTT6l [T31 . In 
these covalently bonded systems at equilibrium, a negative pressure Pbond from the attractive 
covalent bond energy counterbalances the positive one from the repulsion. In Fig. |2l this 
situation corresponds to the 'high pressure' case in which P = |Pbond| an d a negative 
contribution of 5 with v/(2u) < «o < 0. 



Environment-dependent potential for solid Ar 8 

Table 1. Fitted parameters in Vr and p. 

A (J) flitA- 1 ) a 2 (k- 2 ) mi (A- 1 ) M2(A- 2 ) g ^(A" 1 ) 

2.10xl0- 15 -0.5819 0.09309 3.000 -0.03996 80.0 3.60 




Figure 3. Left: Pressure versus lattice constant. Circles denote X-ray observation, (see 
references in [ 1 1) Right: Elastic constants versus pressure: The present model (solid lines) 
compared with experimental results by Shimizu et al. U. 



Finally, parameterisations for functions Vr(R) in Eq.(l20l) and p(R) in Eq.(fT6l) must be 
determined. They are basically similar and described by a superposition of the square of 
two-center overlap integrals. Using Eqs. (fTTT) and (PT41) as a guide, we employ the following 
functions, namely 

V R (R) = A(l + a±R + a 2 R 2 ) exp(-HiR - p 2 R 2 ) (34) 

and 

p(R)=gex V (-uR), (35) 

where the parameters A, p,i,g,v are essential, and ai, a 2 , p>2 are for flexibility of the fitting. 

We do not explicitly include the pairwise E Tep in the present model, because it is actually 
unknown, but may not be dominant, and therefore we may think of it as being absorbed 
effectively in the pairwise component of overlap repulsion unless it proves significant. 



4. Results 

Using the model of overlap repulsion plus the pairwise van der Waals potential (— cqR^ 6 ) 
with the Lennard- Jones parameters for argon ll27ll . i.e. cq = Aea 6 with e = 1.67 x 1CT 21 J 
and a = 3.40 A, the parameters are fitted to the results[8] by ab initio full-potential LMTO 
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calculations with the generalized gradient approximation of PW91 (GGA-PW91)[28] for the 
exchange-correlation, since GGA-PW91 results are remarkably in good agreement with the 
experimental results for solid argon. However, GGA-PW91 is known to predict always 
positive pressures, and a small positive pressure even at the experimental lattice constant 
a = 5.13 Afor the equilibrium at zero pressure. The adjusted parameters are listed in Tabled] 

Figure |3] shows the elastic properties of fee solid argon predicted by the present model 
compared with the experimental results. The agreement is impressive. However, a negative 
curvature in predicted S at very high pressures can be seen as a small deviation, which is 
reflected in B, Cu, C 12 (and not in C 44 ) as it can be understood from 5-term (and its absence) 
in Eqs. ((2~3l) -(l2~6l). This will be corrected if we design more flexible function for p. But if we do 
so, we also need to evaluate neglected terms, for example, the higher order many-atom effects 
due to (1 + 0) _1 factor, which will be handled in a separate study. 

The stability problem of fee against hep is still very subtle even with the present 
environment-dependent model without zero point energy. The result is very sensitive to the 
cutoff. For example, using only the present repulsive model (without E v aw), the fec-hep 
difference in enthalpy is evaluated to be —0.01 meV at 20 GPa and 0.11 meV at 60 GPa if 
we include 86 neighbours within 6 shells in fee and equivalently within 9 shells in hep. The 
difference in the zero point energy [0 would be dominant. Therefore the environmental or 
many-atom effect in overlap repulsion may not be a remedy for the problem of fee stability. 

5. Conclusion 

We have developed an analytic and environment-dependent interatomic potential for the 
overlap repulsion in solid argon. The functional form, of environment-dependence in 
particular, is simple and physically transparent, being based on the non-orthogonal Tight- 
Binding theory for the closed-shell systems. The present model was shown to well reproduce 
the observed elastic properties of solid argon including the Cauchy violation at high pressures. 

A useful and novel analysis has clearly demonstrated how the elastic properties are 
related to the environment-dependence incorporated into a generic pairwise potential. It 
is speculated that the present functional provides not only excellent description for elastic 
properties of a solid noble gas, but also useful description for the problem of small or negative 
Cauchy pressures in covalently bonded systems. 
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